suppressPackageStartupMessages({
  library(readr)
  library(tidyverse)
  library(reshape)
  library(popbio)
})

Qualitative description of population structure

Our data is composed by length data, which we transformed into expected age using the transformed Von Bertalanffy function \[ A = \frac {1}{-K} \times log(1-(\frac {L}{L_{inf}}))+t_0 \], where: A = age, K = growth rate, L = length, \(L_{inf}\) = asymptotic length. Therefore, our model use a age structure.

Description of the demographic information

A life-cycle diagram for your species

DiagrammeR::grViz("
      digraph boxes_and_circles{

# Define nodes
node [shape = circle
      penwidth = 2,
      fontsize = 24]

egg; a1; a2; a3; a4; a5; a6; a7; a8; a9; a10; a11; a12; a13

# Add edge statements
#Advaance classes
egg -> a1[label = s0]
a1 -> a2[label = s1]
a2 -> a3[label = s2]
a3 -> a4[label = s3]
a4 -> a5[label = s4]
a5 -> a6[label = s5]
a6 -> a7[label = s6]
a7 -> a8[label = s7]
a8 -> a9[label = s8]
a9 -> a10[label = s9]
a10 -> a11[label = s10]
a11 -> a12[label = s11]
a12 -> a13[label = s12]


# Fecundities
a1 -> egg[label = f1]
a2 -> egg[label = f2]
a3 -> egg[label = f3]
a4 -> egg[label = f4]
a5 -> egg[label = f5]
a6 -> egg[label = f6]
a7 -> egg[label = f7]
a8 -> egg[label = f8]
a9 -> egg[label = f9]
a10 -> egg[label = f10]
a11 -> egg[label = f11]
a12 -> egg[label = f12]
a13 -> egg[label = f13]
}
      ", height = 1000)

Figure 1 - Life-cycle diagram for Skipjack tuna (Katsuwonus pelamis). S indicates survival, and f indicates fecundity.

Will your model have a pre-breeding census or post-breeding census?

The model we have defined axplicitly models eggs as a cohort (age = 0). We will conssider deleting this in the future, and explicitly calculate recruits.

The matrix written out symbolically

#Build a matrix with text
A <- matrix(0, 14, 14) #initial empty matrix with all 0

# By using the $$, we can embed latex equations that will then be rendered to loog better
# Populate matrix with mortality
for (i in 2:14){
  A[i,i-1] <- paste0("$s_{", i-1, ",", i-2, "}$")
}

ages <- seq(0:13)-1
A[1,] <- paste0("$f_{", seq(1:14)-1,"}$")

knitr::kable(A, col.names = paste0("$a_{",ages,"}$"), row.names = F, caption = "Table 1- Symbolical population matrix A. Subindices represent their rows and columns, respectively. The inferior diagonal represents survivals, while the first row represents facundities.")
Table 1- Symbolical population matrix A. Subindices represent their rows and columns, respectively. The inferior diagonal represents survivals, while the first row represents facundities.
\(a_{0}\) \(a_{1}\) \(a_{2}\) \(a_{3}\) \(a_{4}\) \(a_{5}\) \(a_{6}\) \(a_{7}\) \(a_{8}\) \(a_{9}\) \(a_{10}\) \(a_{11}\) \(a_{12}\) \(a_{13}\)
\(f_{0}\) \(f_{1}\) \(f_{2}\) \(f_{3}\) \(f_{4}\) \(f_{5}\) \(f_{6}\) \(f_{7}\) \(f_{8}\) \(f_{9}\) \(f_{10}\) \(f_{11}\) \(f_{12}\) \(f_{13}\)
\(s_{1,0}\) 0 0 0 0 0 0 0 0 0 0 0 0 0
0 \(s_{2,1}\) 0 0 0 0 0 0 0 0 0 0 0 0
0 0 \(s_{3,2}\) 0 0 0 0 0 0 0 0 0 0 0
0 0 0 \(s_{4,3}\) 0 0 0 0 0 0 0 0 0 0
0 0 0 0 \(s_{5,4}\) 0 0 0 0 0 0 0 0 0
0 0 0 0 0 \(s_{6,5}\) 0 0 0 0 0 0 0 0
0 0 0 0 0 0 \(s_{7,6}\) 0 0 0 0 0 0 0
0 0 0 0 0 0 0 \(s_{8,7}\) 0 0 0 0 0 0
0 0 0 0 0 0 0 0 \(s_{9,8}\) 0 0 0 0 0
0 0 0 0 0 0 0 0 0 \(s_{10,9}\) 0 0 0 0
0 0 0 0 0 0 0 0 0 0 \(s_{11,10}\) 0 0 0
0 0 0 0 0 0 0 0 0 0 0 \(s_{12,11}\) 0 0
0 0 0 0 0 0 0 0 0 0 0 0 \(s_{13,12}\) 0

The matrix filled in with values

Define demographic parameters

# Define Von Bert parameters
# Garbin & Castello, 2014
# L_inf <- mean(c(80, 80, 87.12, 94, 97.9, 97.3, 112.34, 89.38, 92.5))
# K <- mean(c(0.32, 0.6, 0.22, 0.38, 0.14, 0.25, 0.14, 0.38, 0.16))

# Or from Us and Stiurskjgfas, 1981
L_inf <- 102
K <- 0.55
t_0 <- -0.02

# Just to be clear, we use L-inf  and to from Ushisomething, and K from the review

# Define fecundity parameters
# From fishbase http://www.fishbase.se/Reproduction/FecundityList.php?ID=107&GenusName=Katsuwonus&SpeciesName=pelamis&fc=416&StockCode=121
fec_a <- -1.33354
fec_b <- 3.238
# Define mortality

m <- 0.63
z <- 1.39

Define the functions we will need

# Convert length to age using von bertalanffy model, solving for t
length2age <- function(length, l_inf, K, t_o){
  age <- (1/-K)*(log(1-(length/L_inf))) + t_o
  return(age)
}

# Convert age to length using von bertalanffy model
age2length <- function(age, l_inf, K, t_o){
  length <- l_inf*(1-exp(-K*(age-t_o)))
  return(length)
}

#Convert length to fecundity (number of eggs)
fecundity <- function(length, a, b){
  f <- 10^(a+(b*log10(length*10)))
  return(f)
}

Create the matrix

A <- matrix(0, 14, 14) #initial empty matrix with all 0

# Populate matrix with mortality
for (i in 2:14){
  A[i,i-1] <- exp(-1-z)
}

# Populate matrix with fecundity
ages <- seq(0:13)-1
lengths <- age2length(ages, L_inf, K, t_0)
A[1,] <- fecundity(lengths, fec_a, fec_b)

A[is.na(A)] <- 0
A[2,1] <- 0.0000001
A[1,1] <- 0

colnames(A) <- ages
rownames(A) <- ages

knitr::kable(A, col.names = paste0("$a_{",ages,"}$"), row.names = F, caption = "Table 2- Population matrix A. The inferior diagonal represents survivals, while the first row represents facundities.")
Table 2- Population matrix A. The inferior diagonal represents survivals, while the first row represents facundities.
\(a_{0}\) \(a_{1}\) \(a_{2}\) \(a_{3}\) \(a_{4}\) \(a_{5}\) \(a_{6}\) \(a_{7}\) \(a_{8}\) \(a_{9}\) \(a_{10}\) \(a_{11}\) \(a_{12}\) \(a_{13}\)
0e+00 1.657256e+07 7.026737e+07 1.294406e+08 1.758241e+08 2.072322e+08 2.27012e+08 2.38998e+08 2.461086e+08 2.502768e+08 2.527038e+08 2.541114e+08 2.549259e+08 255396708
1e-07 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0e+00 9.162970e-02 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0e+00 0.000000e+00 9.162970e-02 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0e+00 0.000000e+00 0.000000e+00 9.162970e-02 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0e+00 0.000000e+00 0.000000e+00 0.000000e+00 9.162970e-02 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 9.162970e-02 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 9.16297e-02 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 9.16297e-02 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 9.162970e-02 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 9.162970e-02 0.000000e+00 0.000000e+00 0.000000e+00 0
0e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 9.162970e-02 0.000000e+00 0.000000e+00 0
0e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 9.162970e-02 0.000000e+00 0
0e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 9.162970e-02 0

Load the data

Check population size by cohorts through time

Build a population vector

load("size_dist_1980.RData") #This is the n-vector for 1980, loads it to a vector called n

n <- n %>% {
  .$N}%>% 
  c(rep(0, times = 8))

n_0 <- sum(A[1, 2:14]*n, na.rm = T)

n <- c(n_0, n)

n
##  [1] 2.562588e+19 8.674537e+06 3.136387e+11 2.711014e+10 4.391505e+08
##  [6] 3.880114e+06 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [11] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00

Project the population

project <- popbio::pop.projection(A, n, 15)
pop <- project$stage.vectors %>% 
  as.data.frame() %>% 
  mutate(Age = as.factor(seq(0:13)-1)) %>% 
  gather(Year, N, -Age) %>% 
  mutate(Year = as.numeric(as.character(Year))) %>% 
  select(Year, Age, N)

ggplot(pop, aes(x = Year, y = log(N), color = Age)) +
  geom_line() +
  theme_bw()
Figure 2 - Population size through time, represented by ages.

Figure 2 - Population size through time, represented by ages.

Asymptotic analysis of the matrix

Lambda

\(\lambda =\) 1.4729341

plot(project$pop.changes, type = "b", xlab = "Iterations", ylab = "Lambda")
Figure 3 - Convergence of $\lambda$ through time.

Figure 3 - Convergence of \(\lambda\) through time.

Stabe age structure

pop %>%
  filter(Year == 14) %>% 
  ggplot(aes(x = Age, y = log(N))) +
  geom_bar(stat = "identity") +
  theme_bw()
Figure 4 - Stabe age structure, reached after 15 iterations in the model.

Figure 4 - Stabe age structure, reached after 15 iterations in the model.

Sensitivity

sensitivity(A) %>% 
  image2(box.offset=.1)
Figure 5 - Sensitivity matrix.

Figure 5 - Sensitivity matrix.

Elasticity

elasticity(A) %>% 
  image2(box.offset = 0.1)
Figure 6 - Elasticity matrix.

Figure 6 - Elasticity matrix.

Interpretation of the asymptotic results

The \(\lambda = 1.47\) shows that the population is projected to grow at rate of 47%. The result of the sensitivity analysis suggest that is most effective to intervene on the early ages. Since the species mature early (between ages 1 and 2), there is no necessity of conserving the old/big individuals that have higher fecundity, because the younger ones (in quantity) can provide enough eggs to maintain the population.

Discuss what is missing from the model, and how you might add the missing bits

The model assumes the mortality is constant through ages. Better estimates of mortality/survival for different ages would show the real population structure.

The model assumes that the population is density-independent. However, we believe that fecundity, recruitment and survival of this pelagic species is density dependent. Perhaps a first approach wil be to use a Beverton-Holt approach to calculate recruitment, although this assumes density dependence in all the model.

Currently, all organisms are harvested in the model (excluding eggs). An interesting approach would be to set a size limit that reduces total effort on small sizes to natural effort by excluding fishing effort. This could allow us to evaluate the effect of a minimum catch size, management intervention aligned with the observed elasticities.

Environmental variation is another aspect to consider. It is known that recruitment of fish can be affected by extreme environmental conditions, like El Nino. Further analysis could incorporate the effects of climate variability on recruitment.